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ABSTRACT 

This is the first paper of a series on the methods and results of the Active Galactic 
Nuclei In Cosmological Simulations (AGNICS) project, which incorporates the physics 
of AGN into GallCS, a galaxy formation model that combines large cosmological N- 
body simulations of dark matter hierarchical clustering and a semi-analytic approach 
to the physics of the baryons. The project explores the quasar-galaxy link in a cosmo- 
logical perspective, in response to growing observational evidence for a close relation 
between supermassive black holes (SMBHs) and spheroids. The key problems are the 
quasar fuelling mechanism, the origin of the BH to bulge mass relation, the causal and 
chronological link between BH growth and galaxy formation, the properties of quasar 
hosts and the role of AGN feedback in galaxy formation. 

This first paper has two goals. The first is to describe the general structure and 
assumptions that provide the framework for the AGNICS series. The second is to 
apply AGNICS to studying the joint formation of SMBHs and spheroids in galaxy 
mergers. We investigate under what conditions this scenario can reproduce the local 
distribution of SMBHs in nearby galaxies and the evolution of the quasar population. 

AGNICS contains two star formation modes: a quiescent mode in discs and a star- 
burst mode in proto-spheroids, the latter triggered by mergers and disc instabilities. 
Here we assume that BH growth is linked to the starburst mode. The simplest version 
of this scenario, in which the black hole accretion rate M, and the star formation rate 
in the starburst component M^jburst are simply related by a constant of proportion- 
ality, does not to reproduce the cosmic evolution of the quasar population. A model 
in which M, oc Pburst-^'^*burst, where pburst is the density of the gas in the starburst 
and C, — 0.5, can explain the evolution of the quasar luminosity function in B-band 
and X-rays (taking into account the presence of obscured AGN inferred from X-ray 
studies). The scatter and the tilt that this model introduces in the BH-to-bulge mass 
relation are within the observational constraints. The model predicts that the quasar 
contribution increases with the total bolometric luminosity and that, for a given bulge 
mass, the most massive black holes are in the bulges with the oldest stars. 

Key words: galaxies; formation, active - quasars: general 



1 INTRODUCTION 

1.1 The connection between quasars and galaxies: 
observational evidence and open problems 

Several observational and theoretical arguments suggest a 
direct link between quasars and galaxy formation. 



In 'Galactic Nuclei as Collapsed Old Quasars', 
iLvnden-Beiil (^6^ proposed that a quasar phase is part 
of normal galaxy evolution and predicted the presence 
of supermassive black holes in galaxy cores. The detec- 
tion of massive dark objects in the nuclei of nearby 
galaxies has transformed this speculation into an obser- 
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vational fact. The relation between the black hole mass 
and the mass, luminosity and veloci t y dispersion of the 
host s pheroid jM^ggrria^gt al _^2iJE]S££S£g^^^4££— ^ 



2OOOI: iGebhardt et alJ 120001: iMerritt fc Ferraresd l200lt 



Tremain e et alJl2002t iMarconi fc Huntll2003l: iHaring fc Rb3 

2004) together with the chronological agreement between 
the quasar ep och (2^2 — 3) and the ste llar ages of early- 
type galaxies dCattaneo fc Bernardill2'o03l) prove that black 
hole growth mechanisms are directly related to the origin of 
the bulge component. 

This relation opens two kinds of questions. The first one 
concerns the triggering mechanism that activates black hole 
growth (what starts the accretion). The second concerns the 
mechanism that determines the final black hole mass (e.g. 
what terminates the accretion). 

The first attempts to explain the link between quasars 
and galaxy for mation were in the context of the monolithic 
collapse model llEggen et al.lll962D . Matter with low angu- 
lar momentum falls to the centre first forming the galactic 
nucleus and the bulge, while matter with high angular mo- 
me ntum slowly settles i n to a d isc . 

iToomre fc Toornr3 il972h were the pioneers of com- 
puter simulations of galaxy interactions. They argued that 
galaxy mergers can drive a sudden supply of gas into 
the nuclear region while producing a morphological trans- 
formation of spiral galaxies into ellipticals. Smoothed- 
particle-hydrodynamics s imulations of galaxy merge r s con- 
firmed their intu ition (Barnes fc Hcrnauist Il996l Il99lt 
iMihos fc Hernauist 1994. 1996: Springel 2000). The impres- 
sive results of hydrodynamic simulations together with the 
widespread and successful use of the merger model in semi- 
analytic models of hierarchical galaxy formation (see below) 
contributed to a paradigm shift from the monolithic to the 
merger scenario. 

From the point of view of AGN fuelling, the fundamen- 
tal points of the merging scenario are that AGN are fuelled 
with cold gas from the discs of the merging galaxies and that 
black hole growth is conditional to a triggering process. A 
process with a trigger and a short intrinsic duration is con- 
sistent with the episodic nature of AGN (in the monolithic 
collapse model, the brief life time of quasars was attributed 
to the rapid consumption of low angular momentum gas) . 

In relation to the mechanism that determines the black 
hole mass, two scenarios are possible. The first is that black 
hole growth is determined by fuel availability. In this case, 
the quasar switches off when stars have consumed all the 
gas (e.g. Kauffmann_& Hachnclt 2000). In the second sce- 
nario, black ho le growth is self-reg ulated through mechanical 
(jets or winds. lOmma et alJ l200 4| and references therein) or 
radiative (Compton heating. ICiotti fc Ostrikeilll997l 1200 if) 
feedback. In this second case, supermassive black holes may 
be an essential ingr edient of galaxy formation and not just 
a mere bv-product ijsilk fc Reeslliggi : iGranato etHI 12004 
iDi Matteo et ahlbOOsF ! 



1.2 Cosmological models of the formation of 
quasars and galaxies 

The fuelling mechanism, the timing of the active phase in 
relation to the host galaxy evolution and the mechanism 
that limits the accretion are the three main open problems 
in relation to the link between quasars and galaxies. 



200C 



199£ 



There are two ways of approaching these problems. The 
first is to observe and model in great detail a small sample of 
individual objects. The second is to test if specific assump- 
tions can reproduce the population properties of quasars and 
galaxies in a cosmological volume. Semi-analytic models are 

a p owerful tool for t his se cond t ype of study. 

IWhite fc Ree^ lll97^ and lEfstathiou fc ReesI lll988h 
were the pioneer of galaxy formation and quasar forma- 
tion in a cosmological scenario, in which structures form 
through gravitational instability and hierarchical clustering 
of primordial density fluctuation. Semi-analytic models of 
galaxy formation have given a substantial contribution to- 
wards a more detailed comprehension of galaxy formation 
in a cosmological scenario. The field is now so vast that 
it is no longer possible to cite all authors who contributed 
to its development. The most advanced models are those 
of the research groups in Durha m (e.g. [Col^^l 
and references therein}, Munich fKauffmann et al. 
and references there in; iKauff mann fc Hachnclt 200 0), Paris 
( Hat ton et al .l2003fl a nd San ta Gruz dSomerville fc Primackl 
il999HSomerville et aLlliool . Gonceptually, a semi-analytic 
model is structur ed into two steps. In th e first one, Press- 
Schechter theory ("Press fc Schech terll97l) or N-body simu- 
lations (i.e Kauffmami ct al. 199^ are used to follow gravi- 
tational clustering in the dark matter component. The aim 
of this first step is to construct merger trees for a sample of 
dark matter haloes, which are thought to represent the Uni- 
verse. The second step is to use semi-analytic prescriptions 
to follow the physics of the baryons (cooling, star formation, 
stellar evolution, feedback, chemical enrichment, reprocess- 
ing of light by dust, mergers, bar instabilities, etc.) in the 
merger trees. Models in which merger trees are built by using 
N-body simulations, such as the GalIGS (Galaxies In Gos- 
mological Simulations) model by iHatton et alJ ll2003f) , are 
sometimes called hybrid models to distinguish them from 
pure semi-analytic models, in which merger trees are con- 
structed from Press-Schechter theory. 

Quasar m o dellers have followed a similar path 
jGarlberd Il990l: iHaehneh fc ReesI [l993 iKatz ^~sl\ Il994l: 



Maehoenen et al. Il995l: iBi fc Fane Il99'i| iHaiman fc Loe' 



fttanco ct al. '1999"; 'Nitta 1999; 'Mcnou fc Haiman 
;rcival fc Miller .1999.; Haiman fc Mcnou 200* 



Monaco et al. 


I2OOOI; IValageas fc Schaeffed I2OOOI; iGattanec 


2001, 


2001 


Haiman fc Huil I2OOII iHatziminaoglou et al. 


2001 ; 


Martini fc Weinberg! 1200 ll iNath fc Rovchowdhurv 


2002: Volonteri et al."2003a''bHHatziminaoglou et alJl2003f). 



The assumptions in these papers are different, but most of 
them follow the common st rategy of modelling the q uasar 
luminosity function from the lPress fc SchechterHl97'j) mass 
function of dark matter haloes in combination with various 
assumptions for the probability that a halo of mass Mhaio 
at redshift z contains a quasar with luminosity Lqso- These 
studies do not treat the complicated physics of galaxy for- 
mation and AGN fuelling, but show the fundamental scaling 
relations that a physical model has to satisfy. 

The next logical development is 
ing of the two research programmes 
including quasars into semi-analytic 
galaxy f ormation ("Kauffmann fc Hachnclt" 
Haehnelt fc Kauffmann 2000, 2002 
lOi Matteo et al.l 



l2003f) . From an AGN perspective, this 
development allows a more physical modelling and thus a 



2000 
Enoki ct al 



the merg- 
achieved by 
models o f 
2002 
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more physical understanding, while it reduces the number 
of arbitrary assumptions that the modeller can make. E.g. 
merging rates cannot be changed without changing galaxy 
morphologies, colours, luminosity function, etc. Of course, 
that means that more hypotheses and parameters go into 
the model, but the current data allow to constrain them. 
From a galaxy formation perspective, AGN may provide 
an answer to two long standing problems: the overcooling 
problem in massive haloe s (cD galaxies are too many, too 
massive and too blue; e.g. iHatton et a l. 2003 and references 
therein) and the difficulty of reproducing sub-millimetre 
counts without ending up with ga laxies that are too bright 
at th e present cosmic epoch (e.g. iDevriendt fc Guiderdonil 
I2OOOI) . 

1.3 The AGNICS project 

The Active Galactic Nuclei In Cosmological Simulations 
(AGNICS) project started as an extension of the GallCS 
(Galaxies In Cosmological Simula tions) hybrid galaxy for- 
mation model teatton et al.ll2003l) . AGNICS was developed 
to investigate some aspects of the link between quasars and 
galaxy formation in a broad and coherent context. Combined 
models of quasars and galaxy formation are in a qualitative 
agree ment with the cosmic evolutio n of the quasar popula- 
tion jKauffmann fc Haehnefl l2000t) . but do not reproduce 
the quantitative decrease of the number of bright quasars at 
low redshift. S econdly, most of the AGN population is op- 
ticaUy unseen jUeda et alJ 120031: ISazonov et al.ll2004 . but 
the impact of obscured AGN on models of the evolution of 
the quasar luminosity function has not been explored yet. 
Moreover, an interesting by-product of the AGNICS project 
is the possibility to use AGNICS as a tool for generating 
a virtual sky of quasars and galaxies, which may assist the 
planning of forthcoming observations. The virtual sky may 
be used to train automatic recognition algorithms and test 
sources of bias or incompleteness in observational studies. 

The structure of this first paper is as follows. In Section 
2 we recall the main assumptions and features of the Gal- 
lCS galaxy formation rt iodel used to develo p the AGNICS 
model while we refer to lHatton et al.1 ||200J ) for a more de- 
tailed presentation. We give particular importance to the 
origin of galaxy morphologies because a correct modelling of 
spheroids is essential to derive meaningful results for black 
holes and AGN. In Section 3 we describe the method that 
AGNICS uses to model black hole growth and AGN. There- 
after, we distinguish between two main models: the simplest 
possible model, in which the black hole accretion rate and 
the starburst rate are simply related by a constant of pro- 
portionality (we call this the basic model), and the reference 
model (the model which allows us to find the best agreement 
with the data). In Section 4 we show that the basic model 
does not reproduce the low redshift decrease of the comov- 
ing density of bright quasars, but we also show that we can 
modify the basic model by changing the expression for the 
black hole accretion rate and introducing a factor propor- 
tional to a power of starburst density {M, oc p\^^^^^M ^htiist, 
where M. is the black hole accretion rate, M*burst is the 
star formation rate in the starburst component and pburst 
is the density of the gas in the starburst). For C, — 0.5, the 
modified model (the reference model) sati sfies the observa- 
tional constraints deriving from optical JWolf et alJ 120031: 



ICroom et a'i]l2004l) and X-ray iUeda et alJl2003h observa- 
tions of AGN a nd from the scatter in the black hole-to-bulg e 
mass relation llMarconi fc Hundl20o3 : iHaring fc Rixll2004li . 

In the Conclusion (Section 5), we discuss what we can learn 
from the work presented in this article and which questions 
that are still open. The study of quasar hosts is postponed 
to a future publication of the AGNICS series. 



2 THE GALICS GALAXY FORMATION 
MODEL 

GallCS is a model of hierarchical galaxy formation which 
combines high resolution cosmological simulations to de- 
scribe the dark matter content of the Universe with semi- 
analytic prescriptions to follow the physics of the baryonic 
component. The modules that enter the GallCS model are 
the cosmology (to generate the merger trees), the cooling 
model (to follow cooling of hot gas in dark matter haloes), 
the model for the galaxy internal structure and dynamics, 
the merger model (merging rates and effect of mergers on 
galactic morphologies) , the star formation and stellar evolu- 
tion model (including feedback and metal enrichment) and 
the spectral evolution model (stellar spectra, extinction and 
dust thermal emission). We begin this Section by describ- 
ing the N-body simulation used to generate the merger trees. 
Then we progress into explaining the main physical assump- 
tions that enter each functional unit. 

2.1 The dark matter simulation 

2.1.1 Simulation parameters 

The cosmological N-body simulation used to construct the 
merger trees was carried out by using the parallel tree-code 
developed by Ninin (1999). This simulation was run for a 
flat cold dark matter model with a cosmological constant 
of 0.1^ = 0.667. The simulated volume is a cube of side 
Lbox = 100/i^ooMpc, with ftioo = J/o/lOOkms"^ = 0.667, 
containing 256^ particles of mass 8.3 X lO^M© {Ho is the 
Hubble constant). The smoothing length is of 29.29 kpc. 
The cold dark matter power spectrum was normalised in 
agreement with the present day abundance of rich clus- 
ters ((T8 = 0.88). The simulation produced 100 snapshots 
spaced logarithmically in the expansion factor (1 + 2;)^^ from 
z = 35.59 to z = 0. 

2.1.2 Halo identification 

On each snapshot a friend-of-friend algorithm was run to 
identify virialised groups of more than 20 particles. The min- 
imum mass of a dark matter halo is therefore 1.66 x 10^^ Mq . 
For each halo we compute a set of properties, which include 
the position and the velocity of the centre of mass, the ki- 
netic and the potential energy, the spin parameter and the 
tensor of inertia. 

We fit a triaxial ellipsoid to each dark matter halo iden- 
tified in the snapshot. The semi-axes of the ellipsoid are in 
a ratio determined by the ratio between the three different 
components of the inertia tensor. We shrink the ellipsoid 
until the inner region satisfies the virial theorem. The mat- 
ter within this virial volume determines the virial mass. The 
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virial radius is the radius of a sphere whose volume is equal 
to the virial volume and the virial density is equal to the 
virial mass divided by the virial volume. 

When passing from N-body to semi-analytic simula- 
tions, we idealise dark matter haloes as singular isother- 
mal spheres truncated at the virial radius. The informa- 
tion about each individual halo passed to the semi-analytic 
model is therefore contained in three parameters. They are 
the virial mass, the virial density (the mean density within 
the virial radius) and the spin parameter. 

2.1.3 Merger trees 

The merger trees are computed by linking the haloes iden- 
tified in each snapshot with their progenitors at the previ- 
ous time-step. All predecessors from which a halo has in- 
herited one or more particles are counted as progenitors, 
but only one is the main progenitor, the one that has given 
the largest particle contribution. In the same way one can 
identify a main descendant for each dark matter halo. The 
merging histories that we obtain are thus far more complex 
than those constructed from Press-Schechter theory as they 
include the processes of evaporation and fragmentation of 
dark matter haloes. In such situations we assume that all 
baryons remain in the most massive remnant. 

2.2 From haloes to discs 

2.2.1 Gas cooling and infall 

Newly identified haloes are given a mass of hot gas by using 
a universal baryonic fraction of Slb/^o = 0.135. Hot gas is 
assumed to be shock heated to the virial temperature of the 
halo and in hydrostatic equilibrium with the dark matter 
potential. The density profile of the hot gas is that of a sin- 
gular isothermal sphere truncated at the virial radius. The 
cooling time is calculated from the hot gas density distribu- 
tion with the [Sutherland fc Dopital 1^993) cooling function. 
The gas that can cool in a time At is the gas whose cooling 
time and free-fall time are both lower than At. GallCS as- 
sumes that gas cooling is accompanied by a simultaneous gas 
inflow in order to maintain the same hot gas density profile 
at all times. The disc radius is = Arvir/%/2, where r^ir is 
the virial radius and A is the angular momentum parameter 
of the halo. Cooling is inhibited in haloes with A > 0.5. 

We also prevent cooling in haloes that contain a total 
mass of bulge stars > 10^^ Mq not to exceed the observa- 
tional constraints on the number of massive galaxies. The 
theoretical argument behind this solution is that the bulge 
mass is proportional to the black hole mass and therefore 
to the total energy budget of AGN over the life of the halo. 
We assume that AGN deposit energy in the intergalactic 
medium mechanically (through jets or winds) and radia- 
tively (through hard X-ray photons) and that this input 
is proportional to the black hole accretion rate. In reality, 
the long term outcome of the interaction between the black 
hole and the intergalactic medium is poorly understood. In 
particular, there is no strong physical reason why feedback 
should become important above a critical black hole mass 
(but see .Cattan eo 2Q02; Dunlon ct al. 2003') . The cut-off at a 
bulge mass of > lO" Mq introduced bv lHatton et al.l 1)200,^ ^ 
and used for the work presented in this article is a temporary 



measure, whose justification is in its capacity of reproduc- 
ing the data. In fact, all groups working on semi-analytic 
models of galaxy formation were forced to make arbitrary 
assumptions to fit the bright part of the galaxy luminosity 
function. Meanwhile, we are looking for a motivated self- 
consistent approach, which will appear in another paper of 
the AGNICS series. 

2.2.2 Star formation and stellar evolution 
The star formation rate in the disc is 

M. = (1) 

Here Mcoia is the mass of the gas in the disc (all the gas in 
the disc is cold and all the gas in the halo is hot) and tdyn is 
the dynamical time (the time to complete a half rotation at 
the disc half mass radius). The parameter /3, , which deter- 
mines th e efficiency of star form ation has a fiducial value of 
/3* = 50 llGuiderdoni et al .1199811 . The mass of newly formed 
stars is distributed according to the lKennicutd lll983ll initial 
mass function. Stars are evolved between time-steps using a 
sub-stepping of at most 1 Myr. During each sub-step, stars 
release mass and energy into the interstellar medium. Most 
of the mass comes from the red giant and the asymptotic 
giant branch of stellar evolution, while most of the energy 
comes from shocks due to supernova explosions. In GallCS, 
the enriched material released in the late stages of stellar 
evolution is mixed to the cold phase, while the energy re- 
leased from supernovae is used to reheat the cold gEis and 
return it to the hot phase in halo. The reheated gas can also 
be ejected from the halo if the potential is shallow enough. 
The rate of mass loss in the supernova-driven wind that flows 
out of the disc is directly proportional to the supernova rate. 
The details of the m odel and their justiflcation are found in 
iHatton et aP J2003l) . 

2.3 Galaxy morphologies 

2.3.1 Galaxy internal structure 

GallCS models the baryonic part of a galaxy as the sum 
of three components: the disc, the bulge and the starburst. 
These three components are not always present at the same 
time. The fundamental assumption is that all galaxies are 
born as discs at the centre of a dark matter halo. The trans- 
formation of disc stars into bulge stars and of disc gas into 
starbursting gas is is due to bar instabilities and mergers. 
Gas is never added to bulges directly and the only gas in 
bulges is that coming from stellar mass loss. The starburst- 
ing gas forms a young stellar population that becomes part 
of the bulge stellar population when the stars have reached 
an age of 100 Myr. We do not readjust the bulge radius when 
this happens. The disc has an exponential profile, while the 
bulge and the starburst are described by ajHernquist ( 199^ 
density distribution. The starburst scale is rburst = Krbuigc 
with K — 0.1. This model cannot take into account the pres- 
ence of more complex morphologies, which are particularly 
relevant to the case of starburst galaxies, but has a justifica- 
tion that comes from smoothed particle hydrodynamics sim- 
ulations of galaxy mergers. The star formation law (Eq. 
has the same form and uses the same efficiency parameter 
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/3, for all three components when we redefine M^oid as the 
mass of the geis in the component and tdyn as the dynamical 
time of the component. For the components described by a 
Hernquist profile, the dynamical time is tdyn = '"0.5/0", where 
ro.5 is the half mass radius and a is the velocity dispersion 
at the half mass radius. 

2.3.2 Disc instabilities 

The formation of a bulge from a disc instability is mod- 
elled as follows. A disc is globally stable if Vd < O.Twtot (e.g. 
[van den Bosch 1998). Here wtot is the circular velocity at the 
disc half mass radius, computed from the total gravitational 
potential of the disc, the bulge, the starburst and the halo, 
while Vd is the circular velocity computed from the gravita- 
tional potential of the disc only. The bar instability transfers 
gas from the disc to the starburst and stars from the disc to 
the bulge until the stability criterion is fulfilled. 



gas is transferred to the starburst and stars are transferred 
to the bulge. The fraction of the disc mass transferred to 
the spheroidal component (the bulge and the starburst) de- 
pends on the mass ratio of th e merging galaxies. T he exact 
formula for this fraction is in lHatton et alJ ll2003h . but the 
idea is that when one of the two galaxies is much less massive 
than the other one (minor merger), most of the disc mass 
stays in the disc. On the contrary, if the two galaxies have 
comparable masses (major merger), most of the disc mass 
is transferred to the spheroidal component. The separation 
between a minor and a major merger is for a mass ratio of 
1:3. 

Devriendt et al. (in preparation) have used GallCS to 
investigate the relative importance of bars, minor mergers 
and major mergers in the formation of spheroids. Ellipti- 
cal galaxies acquire most of their mass through mergers, 
although disc instabilities are predominant in forming the 
bulges of spiral and lenticular galaxies. Minor mergers pre- 
dominate over major mergers in all galaxy types. 



2.3.3 Merging rate 

Mergers are the second path to form spheroids. The model 
for galaxy mergers is articulated into two parts: the compu- 
tation of the merging rate and the model for predicting the 
outcome of mergers when they occur. Let us begin with the 
first one. In the beginning there is one galaxy at the centre 
of each halo. When two haloes merge, the central galax- 
ies of the two haloes will sink towards the common centre 
of mass owing to dynamical friction. That will take some 
time, so that the halo may temporarily be without a cen- 
tral galaxy. The common centre of mass will be closer to the 
centre of mass of the more massive halo, so that the cen- 
tral galaxy of the more massive halo will reach the centre 
of the new halo more rapidly, not only because it is more 
massive and is therefore subject to stronger dynamical fric- 
tion (computed from the standard differential equation in 
iBinnev fc Tremaind[l987l) . but also because it starts from 
a lower radial coordinate. Galaxy groups and clusters form 
because above a certain mass the dynamical friction time- 
scale for galaxy mergers becomes longer than the time-scale 
on which dark matter haloes merge. When groups or clusters 
merge, the radial coordinates of satellite galaxies in the new 
halo are determined according to a prescription such that, 
if the satellites come from the main progenitor and this is 
almost as massive as the final halo, their radial coordinates 
are nearly unperturbed. On the other, if the satellites come 
from a much less massive halo, their final radial coordinates 
are close to the virial radius of the new halo. GallCS con- 
siders not only mergers at the centre of the halo due to 
dynamical friction, but also mergers due to encounters of 
galaxies at a non zero radial coordinate (collisions between 
satellite galaxies), but these events are a small correction to 
the predominant merging activity with the central galaxy. 



2.3.4 Effect of mergers 

We now come to the second part of the merging model: 
how mergers transform galactic morphologies. Bulge mate- 
rial remains bulge material and starburst material remains 
starbust material, while some of the disc material is re- 
moved from the disc. As in the case of the bar instability. 



2.3.5 Tests of the model 

In Fig. 1, the predictions of the GallCS model for the mass 
function and the stellar ages of local spheroids are shown as 
solid lines. These predictions are compared with a variety of 
data (shown as points with error bars, dashed and dotted 
lines) . 

The mass function of spheroids is a key test of the 
model. We have tried to make the comparison as robust 
as possible by using data derived with completely different 
methods. The first is to start from the galaxy luminosity 
function decomposed by morphological type, to assume a 
bulge fraction for each morphological type and finally to 
use a bulge mass-to-light ratio to convert the spheroid lu- 
minosity function into a mass function. We have applied 
this method to the 2Mass K band data of iKochanek et aH 
(.2001) and to the Sloan Digital Sky Survey r* band data 
of lNakamura et a k| l|200 2|). both with an aperture correction 
of -0.2 mag fsee lsi^mkar et alJl2004l from where we have 
taken the bulge-to-total mass ratios for early-type and late- 
type galaxies, and the K band and R band mass-to-light 
ratios). With these data, we ha ve found the mass func tions 
shown by the o pen diamonds f Cochanek et al.l [20011) and 
the dotted line IINakamura et al ]200^r Tn Fig. 1. The sec- 
ond way to compute the spheroid mass function is to per- 
form a bulge to disc decomposition on a complete galaxy 
sample, to determine a bulge mass function and then to 
use a bulge mass-to-light ratio. The open triangles show 
the result of this method with the / band spheroid lumi- 
nosity function of Benson et al. (2002). The third method 
does not pass through galaxy luminosities at all. Instead, 
it uses the relation of bulge mass and velocity dispersion 
to convert the Sloan early-type galaxies velocity disper- 
sion distribution (e.g. ICattaneo fc Bernardil 120031) into a 
spheroid mass function. The spheroid mass function com- 
puted through this method is shown by the dashed line 
(which misses the bulges of spirals because it is derived 
from a s ample of early- type ga laxie s). The luminosity func - 
tions bv lKochanek et all i200lh and Nakamura e t al.l ll2002h 
give very simi l ar res ults. The mass functions obtained from 
iBenson et all (|20o3) and ICattaneo fc Bernardil ll2003l) re- 
semble each other very closely, too, while there is difference 
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Figure 1. i He properties ot spJieroids m tnc UaliUb model and m Vac data. 1 He lett pa , 
by GallCS (solid line) with those inferred from the 2Mass K luminosit y function of | Kochanek et alJ l2^^1^ (diamonds with error bars), 
the Sloan Digital S ky Suryey r* band luminosity function of Nakamu ra et all 12003) (dotted line), the I band l uminosity function of 
iBcnson e t al.l l2002^ (triangles with error bars) and the Sloan Digital Sky Suryey velocity dispersion distribution in lCattaneo fc Bernardil 
(2003) (dashed line). The right panel shows the mass - stellar age relation in GallCS (small points; the solid line shows the mean yalue) 
and in the Sloan Digital Sky Survey iCattaneo fc BernardilEoOSl points with error bars). 



between the first and the second group. However, the over- 
all impression is that there is good agreement between the 
model and the data. In this comparison, we should have in 
mind that the dark matter simulation cannot resolve haloes 
less massive than 1.6 x lO^^Af©. (haloes that contain less 
than 2 x lO^^M© in baryons). Therefore, a galaxy less mas- 
sive than 2 x 10^" M0 is formally below the resolution limit. 
In the case of Fig. 1, the importance of resolution effects is 
reduced by the fact that the fraction of the baryonic mass 
in bulges is small in low mass haloes. 

The second panel of Fig. 1 shows the mean stellar ages 
of spheroids in GallCS (solid line) and in the Sloan Digi- 
tal Sky Survey (points with error bars). The small points 
show the scatt er in GallCS age s. The points with error 
bars come from lCattaneo fc Bern ardi (2003) using Afbuigc ~ 
10^^ (cr/200 kms~^)*, where a is the stellar velocity disper- 
sion of the spheroid. The stellar ages of spheroids that we 
are comparing are not exactly the same thing because the 
model outputs are mean values weighted by stellar mass, 
while the mean values inferred from the Sloan Digital Sky 
Survey are weighted by luminosity, but the difference should 
not be large for the r-band, especially for massive galaxies 
with old stellar populations. The agreement is acceptable 
but not so good as that found for the spheroid mass func- 
tion. In the cosmology used for this paper the Universe is 
13.7 Gyr old. Lookback times of 7, 10 and 11 Gyr correspond 
to redshifts of 1, 1.7 and 2.3, respectively. Low mass ellipti- 
cals are older in the model than in the data because at low 
masses disc galaxies contaminate the early-type sample, but 
the real problem is at large masses. GallCS shows that the 
typical redshift of formation of the stars in a galaxy with 
Mbuigc — 2.5 X lO^^M© is of 2 ~ 1.1 whereas the data sug- 



gest a value closer to z — 1.7. This discrepancy is a sign 
that the formation of spheroids in semi-analytic models of 
galaxy formation is an open problem and indicates the need 
of more work on this topic. 

In Fig. 2 we show the cosmic evolution of the galaxy 
merging rate by plotting total masses (stars and cold gas) 
and gas masses for all mergers that have occurred in 4 time 
spans of 100 Myr each at redshifts of about 2, 1, 0.5 and 
0. This figure contains more information than a plot of the 
merging rate as a function of redshift for galaxies above a 
mass threshold. The decrease of the merging rate as a func- 
tion of redshift is not dramatic because the larger intergalac- 
tic distance at lower density of is compensated by a larger 
galaxy number, especially at high galaxy masses. On the 
other hand, low redshift mergers have a tendency to a lower 
gas fraction. 

2.4 Galaxy luminosities 

Stellar spectral energy distributions (SEDs) in the optical 
and near infrared are computed by convolving the star for- 
mation history of each component with the S EDs derived 
from the STARDUST jPevriendt et al.lll999^ stellar pop- 
ulation synthesis model. STARDUST calculates the SEDs 
associated with a single burst of star formation at time in- 
tervals that go from 10 Myr to 50 Gyr. These SEDs assume 
a Kennicutt stellar initial mass function and depend on the 
metallicity of the stellar population. Dust absorption is com- 
puted with an extinction law depending on the mean hy- 
drogen column density and the gas metallicity. The column 
density depends on the mass and the geometry of the gas 
distribution. GallCS uses two models: a spheroidal distribu- 
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Figure 2. Total masses (stars and cold gas, x-axis) and gas masses only (t/-axis) for mergers that have occurred in four time intervals 
of 100 Myr each at redshifts of 2, 1, 0.5 and 0. The dashed lines correspond to the position on the plots of mergers in which the galaxies 
are 100% gas. Notice that the plotting range on the vertical axis is different in each plot. 



tion for bulges and starbursts and a uniform slab for discs. 
In both cases, stars and dust have the same space distribu- 
tion. For each disc, GallCS picks a random inclination angle 
and computes the extinction for that value. All the radiation 
that is absorbed is re-emitted in the infrared by four dust 
components: big carbonaceous grains, small carbonaceous 
grains, silicates and polycyclic aromatic hydrocarbons. The 
infrared colour-luminosity relation observed in IRAS galax- 
ies determines how the absorbed power is distributed among 
the four components. 



3 THE STRUCTURE OF THE AGNICS 
MODEL 

3.1 Black hole growth model 

The most recent data suggest that most bulges contain a 
supermassive black hole (see the Introduction) . We do not 
know how these black holes formed, although iReed (Il984l) 
identified a number of astrophysical paths that are likely 
to result in the formation of a supermassive black hole. 
The main problem is how gas concentrates from a scale 
of ~ 100 pc down to ~ 1 pc. We do not try to model the 
complicated physics of how gas flows from the galaxy into 
the accretion disc and then from the accretion disc into the 



supermassive black hole. Instead, we accept as a fact that 
whenever there is a bulge or a starburst, there is a super- 
massive black hole. 

The model for black hole growth is based on three in- 
gredients: the initial black hole mass, the model for black 
hole accretion and the model for black hole coalescence. 

3.1.1 Seed mass 

The initial black hole mass is important in models where the 
black hole mass determines the maximum accretion rate. If 
the accretion rate cannot exceed the Eddington limit (dis- 
cussed in Section 3.2), then a larger initial mass allows the 
black hole to grow more rapidly. If we do not limit the ac- 
cretion rate to the Eddington limit, then the value of the 
initial black hole mass is irrelevant (at least as long as the 
initial black hole mass is M. < 10^ M©). 

3.1.2 Accretion rate 

The model for black hole accretion is the fundamental ele- 
ment that distinguishes between different models in the lit- 
erature and also between different implementation of the 
AGNICS model. There are two families of models: those 
in which black holes grow through the accretion of cold 
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gas (e.g-ICattaneo et ailligogtlKauffmann fc Haehneltl200ol: 
lEnoki et al.ll2003ll . in which case we expect a relation be- 
tween black hole growth and star for mation, and t hose in 
which black hole gr ow through spherical | Bondi lIlQSST) accre- 
tion of hot gas fe.g. lNulsen fc Fabianl200(]ll . In this paper we 
only fuel AGN with cold gas. We tried to develop models of 
the second type but we were not able to do it without over- 
predicting the number density of bright quasars at low red- 
shift. The reason is that, at low redshift, potential wells are 
deeper and haloes are less dense, with the consequence that 
shock heating is more effective, the cooling time is longer 
and a higher baryonic fraction is in hot gas. 

Within the family of models in which black hole ac- 
cretion is fuelled with cold gas and triggered by mergers, 
we can identify different p rescriptions for t he g as mass ac- 
creted by the black hole. ICattaneo et alJ il999,i'l assumed 
that at each major merger the black hole accretes a mass 
AM, oc (1 -I- 2:)''AM,burst, where AM«burst is the stellar 
mass formed by the merger (in practice the mass of the 
gas in the discs of the merging galaxies) and t; is a free 
parameter determined to reproduce the cosmic growth of 
the comoving mass density in supermassive black holes (a 
reasonable agreement with the data required that 77 ~ 2). 
[Kauffmann & Haclmclt (2000) assumed that on average, at 
each major merger, the black hole accretes a mass 

/bhMcom 



AM. 



(2) 



l + (280kms-VT4)2' 

where Mcoid is the mass of cold gas in the merging galaxies, 
Vc is the host halo circular velocity and /bh is a free param- 
eter determined from the black hole-to-bulge mass relation 
The accretion was distributed over time with the law: 



AM. 

M. = exp 

iaccr 



V iaccr 



(3) 



where taccr(2) is a free para meter determined b y fitting the 
quasar luminosity function. lEnoki et alJ (|200^ used a very 
similar model with 

_ /bhMcom 



AM. = /bhAA/.; 



(4) 



0.75 + (280kms-VK)2-^ 

where AM,burst is the total mass of the stars formed in the 
starburst that accompanies the merger. 

In AGNICS we integrate a differential equation for the 
black hole mass. In this paper we only consider models of 
the form: 



M. = fcM.burst. 



(5) 



M*burst is the star formation rate in the starburst compo- 
nent. Black holes in starbursts are active while black holes in 
galaxies that do not contain a starburst are quiescent. The 
accretion efficiency k can depend on an large number of pa- 
rameters (morphologies of the merging galaxies, gas fraction, 
structure of the intergalactic medium, effect of AGN and su- 
pernova feedback, impact parameter, inclination of the or- 
bital plane with respect to the discs of the merging galax;ies, 
corotation or counterrotation of the two discs, black hole 
spin, etc.). In practice, we only consider very simple cases. 
In this paper we discuss four of them (Table 1 summarizes 
all the models considered in this paper): 

• k =constant: there is a one-to-one relation between 
black hole growth and starbursts. This is what we called 
the basic model and corresponds to model A. 



• k ^constant, Af. < MEdd (model B). 

• fc oc /5^oi;j where pcoid is the density of the gas in the star- 
burst component. This is the reference model, which gives 
the best fit to the data and corresponds to models C, D and 
E. 

• fc oc -f (280kms^"'^/Vc)^] to suppress black hole 
growth in small haloes (model F). 

These recipes describe a mean accretion rate: the ap- 
proach model oversimplifies the physics of individual ob- 
jects, but can be expected to be meaningful in a statistical 
sense, when it is used to compute mean values in a cos- 
mological volume. However, neglecting the dependence of k 
on other parameters will lead to underestimate the scatter 
in the properties of black holes and AGN. The scatter in 
the results of the simulations can only be meaningful as a 
lower limit to the scatter predicted by the model. Another 
remark: Eq. implies that the AGN duty cycle coincides 
with the entire duration of the starburst. This is a simpli- 
fication because black hole accretion is likely to be highly 
time- dependent, as it is shown by AGN variability. 



3.1.3 Black hole coalescence 

The third choice that one has to make to complete the black 
hole growth model is what happens when two galaxies with 
black holes merge. We assume that the black holes imme- 
diately coalesce prior to any gas accretion. This is not a 
good description of reality if the active phase starts before 
the black holes have merged (it is not difficult to find im- 
ages of quasar hosts with double nuclei). It is also possi- 
ble that a third black hole from a second merger reaches 
the centre before the two black holes from the first merger 
have had time to coalesce. In this case, the three-body 
interaction may result in the ejection of one black hole 
from t he galaxy teegelman et alJIlOSd : iHut fc Ree'3ll992l: 
IVolont cri ct al. 2003b). Nevertheless, the simple assumption 
of instantaneous coalescence is the most natural within the 
assumption of instantaneous morphological transformation 
at the time of merging used in semi-analytic models. 



3.2 Luminosity and spectral energy distribution 
of AGN 



The bolometric luminosity of a quasar is 



^rad 



^rad 



-M.c" 



(6) 



where e^ad ~ 0.1 is the radiative efficiency of black hole 
accretion. Bardccn ( 1991) explains the physics of this value. 

The Eddington luminosity Lsdd is the critical lumi- 
nosity above which the radiation pressure is stronger than 
the gravitational attraction. The justification for requiring 
Lhoi < -f'Edd is that, as soon as the Eddington luminos- 
ity is exceeded, radiation pressure pushes the gas outwards 
and brings the accretion rate back to sub-Eddington val- 
ues. However, there are physical situations in which this 
limit can be exceeded, e.g. if the emission is beamed into 
a narrow cone or if photons are trapped into a radiatively 
inefficien t, optically thick flow and advected into the event 
horizon I^Begelmanlll978^ . In the latter case it is possible 
to have M, > A^Edd without requiring that Lboi > ^Edd 
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(here MEdd is the black hole accretion rate for which Eq.|S| 
gives I/boi = L-Edd)- Imposing the condition M, < M^dd 
introduces a characteristic time-scale for the growth of the 
black hole. This time-scale is linked to the Salpeter time 
ts = 4 X 10* yr, the time in which a black hole radiating at 
the Eddington limit releases its entire mass energy M,c^. 
Since a black hole cannot radiate its entire mass energy, the 
accretion time-scale is limited to taccr = tradts and the black 
hole mass cannot grow faster th an exp {—t/ tg^c.T ) • 

'Mi rabel fc Rodrfgued (|l99^ have shown that some mi- 
croquasars in the Milky Way are accreting well above the 
Eddington limit. In this paper we consider different models 
(see Table 1): the default assumption is that nothing stops 
black holes from accreting at super-Eddington rates, but the 
luminosity cannot e xceed the Eddington limit. Th is is the 
same assumption of lKauffmann fc Haehnell ll200(J) . 

The sta ndard reference in the literature for the quasar 
SED is the lElvis et alJ median SED inferred from 

radio, sub-millimetre, infrared, optical, ultraviolet and X- 
ray observations of a sample of optically selected quasars. 
With this template, the blue magnitude of an optical quasar 
is 

Mb = -23.64 - 2.5 Log(Lboi/10*'' erg s"') (7) 

and the colours are: U - B = -0.9 11, B-V = 0-092 , V - 
R = 0.320 a.ndR-K ^ 3.287. The lElvis et"ai] (|l993) SED 

gives a bolometric co r rectio n of ~ 32 for the 2-10 keV X-ray 
band. iMarconi et al.l ||2004') have shown that a luminosity- 
dependent bolometric correction provides a more adequate 
fit to the X-ray data and it is their model that we use to 
compute X-ray luminosities. 




^ 0.2 - .... 

n f) r I I I .. I ... .1. 1 ' 

41 42 43 44 45 46 
Lx (erg/s) 2-10 keV 

Figure 3. The fraction of X-ray selected AGN that do not have an 
AGN op tical/UV continuum or broad emission lines. The squares 
are fro m lUeda et al.l 1200,'^) . the t riangles are fr om Szokolv et all 
1*200?) and the diamonds are from lBarger et al.l (2005'l. The solid 
line (extinction law I) and the dashed line (extinction law II) are 
two models used to fit these data. 



3.3 Obscured AGN 

Obscured AGN are AGN where the optical/UV continuum 
and the broad line spectrum are not observable. Unified 
models state that unobscured and obscured AGN are in- 
trinsically the same type of systems, made of a central black 
hole surrounded by its accretion disc and, more further out, 
by a dusty torus. In obscured (also called type 2) AGN, the 
system is observed from the equatorial plane of the torus and 
the optical/UV continuum is not observable because it is ab- 
sorbed by dust on the line of sight. Obscured AGN can be de- 
tected in hard X-rays, at optical wavelengths (through scat- 
tered polarised light and strong narrow line emission), by 
using mid-infrared spectroscopy and with radio observations 
(see the review article by Ris a liti fc Elvis 20q4').IUeda et al.l 
l2003h : ISzokolvet al.l ^l2004^ : iBarger et alJ ^2005^ used X- 
ray selected sample to estimate the fraction of obscured 
AGN and found that the result is strongly luminosity- 
dependent (Fig. 3). Most bright quasars are unobscured, 
while > 75% low luminosity AGN are obscured . These 
findings are consistent with those bv lSazonov et all i2004l) . 
who studied the combined cosmic infrared and X-ray back- 
ground and inferred an obscured/unobscured ratio of a fac- 
tor of three. Optical extinction is therefore an important 
effect and can shift the peak of the cosmic accretion his- 
tory of supermassive black holes to lower r edshifts (e.g. 
ICatt anco fc Bcrnardi '2003': "Stc ffen et a'i] |2003^ . 

The Elvis ct al. ( 1994) SED is a phenomenological tem- 
plate for unobscured (type 1) quasars. It contains three com- 
ponents: the big blue bump (due to thermal emission from 



the accretion disc), the infrared bump (due to thermal emis- 
sion from the inner edge of the dusty torus) and the X-ray 
bump (due to Comptonization of accretion disc light). We 
cannot treat AGN extinction as we do with the absorption 
of stellar light by dust, where we assume a covering factor of 
~ 1, because that grossly underpredicts the number of AGN 
by extinguishing the entire quasar population with Ay ~ 5. 
Instead, we assume that for a given bolometric luminosity 
there is an obscuration probability fohs{Lx), which is de- 
termined phenomenologically from the data in Fig. 3. This 
is the same as to assume that the torus has a luminosity- 
dependent opening angle. When an AGN is viewed from the 
equatorial plane of the torus we still see the infrared ther- 
mal emission from the torus and the X-ray emission from the 
central engine (assuming that we can neglect self-absorption 
at infrared wavelengths and that the torus is not Compton- 
thick). 

The most delicate point is the choice of the function 
fobs- Fig. 3 shows the two model functions used in this 
pape r. The first one, show n by the dashed lin e, fits the 
data lUeda et alJ lj2003h and lSzokolv et al.l (|2004h . while the 
second o ne, shown by the s olid line, is closer to the new 
data by iBarger et alJ JioO^. The two lines cover only a 
small fraction of the luminosity range because AGN with 
Lx <J 3 X 10*^ ergs~^ are too faint to contribute to the lu- 
minosity function of optical quasars and we do not want to 
extrapolate the curves beyond the range to which we are ap- 
plying them as a phenomenological model for the extinction 
of quasars by dust. 
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Model k Eddington limit e^ad Extinction 

A (the basic model) 0.0017 
B 0.0017 



C (the reference model) 0.0012*^0^5 r"!^'^ 

D 0.0012M0-5rr^"^ 
E 0.0012M0-5r"^'5, 

11 burst 

F 0.0012M0i5r.-^l^'5/[l + (280kms-l/K)2] 



No 


0.26 


II 


M, < J\/Edd 






Lbol < -f-Edd 


0.1 


II 


^bol < ^Edd 


0.1 


I 


No 


0.08 


II 


Lhol < -f-Edd 


0.1 


I 



Table 1. Models used to produce the figures of the paper. In the most general case, the black hole accretion efficiency k (Eq.|2l depends 
on the mass of the gas in the starburst Mn (in units of 10^^ Mq), on the starburst radius rtm-st (in kpc) and on the halo circular 
velocity V28O (in units of 280kms~^). We have considered three models for the physics of the Eddington limit. In models A and E, the 
Eddington limit is completely ignored. In model B, black holes are constrained to accrete at sub-Eddington rates. In models C, D and 
F, black holes can grow at super-Eddington rates but their luminosities cannot exceed the Eddington limit. The radiative efficiency e,.^^ 
is the one used in Eq. JHJ to convert the black hole accretion rate into a bolometric luminosity. The extinction laws I and II refer to the 
solid line and the dash line in Fig. 3, respectively. 



Model A z = 




Model B z = 




10.0 10.5 11.0 11.5 12.0 
Log(M,„,ge/M0) 



10.0 10.5 11.0 11.5 12.0 
Log(M,„„e/Mo) 



Figure 4. The simulated black hole mass - bulge mass relation in a model in which the black hole accretion rate is directly proportional 
to the star formation rate in the starburst component (M, = O.OOlTMtburst)- The left and the right panel shows the model predictions 
respectively without and with the condition M, < A/gdd (models A and B of Table 1). Model B is computed for a seed mass of 
M, = lO^'^Mcoid where M^oid is the mass of the gas in the starburst component. The predictions made by AGNICS are shown as dots, 
bl ue for bla ck holes in disc-dominated galaxies, red for black holes in bulge-dominated galaxies. The open squares are the mass estimates 
by iMarcon i & Hunt 1 200^ and the filled triang les those by Haring & Rix 12004). The dashed-dotted lines corresponds to a 1.12 slope 
(th e best fit to the data of H aring fc Rixll2004l) and the dashed lines show a scatter of ±0.3 dex around this relation. The solid line is 
the lMarconi fc Hunh 1.2003.) M, = 0.002Mbuigc relation rescaled for M, = 0.0017M]-,uige. These lines and the data points are the same 
in both panels. 



BLACK HOLE MASSES AND THE 
REDSHIFT EVOLUTION OF THE QUASAR 
POPULATION 



We want to determine if a model in which the black hole 
accretion rate is related to the star formation rate in the 
proto-spheroid (Eq. |3 can reproduce the masses of black 
holes in nearby galaxies and the cosmic evolution of the 
quasar population. 



4.1 The basic model 

The simplest scenario is that the black hole accretion rate 
is directly proportional to the star formation rate in the 
starburst component (models A and B of Table 1). This as- 
sumption gives a tight M, - Mbuigc relation consistent with 
a linear scaling. The left and the right panel in Fig. 4 present 
the case without and with the constraint M, < MEdd, re- 
spectively. The blue and the red points are the results of 
the simulations for disc-dominated (late-type) and bulge- 
dominated (early- type) galaxies. The ope n squares are the 
mass estimates bv lMarconi fc Hunti ||20o3) and the filled tri- 
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angles those bv lHaring fc Rij i2004ll . The sohd hne shows 
the hnear relation M, = 0.0017Mbuigc- These lines and the 
data points are the same in both panels. 

The tight M, oc Mbuigc relation follows directly from 
our assumption. The ~ 0.1 dex scatter in the left panel 
(model A) is due to the presence of a stellar component 
in the discs of the merging galaxies. In a merger, the masses 
of the two bulges add together, and so do the masses of the 
two black holes. The gas in the discs of the merging galaxies 
contributes to the formation of bulge stars and to the growth 
of the black hole in a fixed proportion. However, the stars in 
the discs of the merging galaxies can only contribute to the 
growth of the bulge mass, and as the gas fraction is lower 
at high masses (Fig. 2), high mass mergers produce lower 
values of M,/Mbuigc and a shallower M, - Mbuigo relation 
at high masses. The model with M, < Msdd (model B) con- 
tains more scatter because the Eddington limit introduces 
a characteristic time-scale for the growth of the black hole. 
In some galaxies with a short star formation time-scale, the 
starburst may be over before the black hole has had time to 
grow significantly. Nevertheless, the similarity between the 
two panels of Fig. 4 indicates that the importance of this 
effect is limited. 

The basic model cannot reproduce the observational 
scatter in any of the two versions. That is not surprising, 
because we expect that our approach underestimates the 
scatter, since it ignores many parameters that can affect the 
mass accreted by black holes (Section 3.1). The real problem 
comes from the quasar luminosity function (Fig. 5). 

Once we have specified the black hole growth model, 
our only freedom to fit the quasar luminosity function is in 
the value of the radiative efficiency parameter trad and in the 
choice of the extinction law (the fraction of type 2 AGN as a 
function of the X-ray luminosity, /obs(ix)). The luminosity 
functions corresponding to models A and B in Fig. 5 were 
computed for a radiative efficiency of trad — 0.26 by using 
the extinction law II in Fig. 3. Very likely, these assumptions 
overestimate the blue light that comes out of powerful AGN. 
Nevertheless, Fig. 5 shows that model A underpredicts the 
number of bright quasars at 2; ~ 2 (model B is even worse) . 

We can improve the fit to the luminosity function at 
z ~ 2 by pushing the radiative efficiency to an even higher 
value, by removing optical extinction completely and by in- 
creasing the normalisation of the M, - Mbuige relation (the 
data allow an increase up to a maximum of 30%). How- 
ever, that misses the physical point, besides the difficulty 
of justifying such assumptions. The problem exists because 
the basic model underestimates the scatter in the black hole 
mass distribution. If the black hole mass distribution is a 
Gaussian in Log(M,), then more scatter in Log(M.) with 
the same mean value of Log(A/.) increases the comoving 
mass density of supermassive black holes, which automati- 
cally raises the total light output of the quasar population. 
Therefore, we can fix the lack of bright quasars at z ~ 2 just 
by introducing a scatter factor cx lO'^ in front of the right 
hand side of Eq. (|KJ , where x is a Gaussian random deviate. 

However, that does not solve the second problem of 
model A: the simulated evolution of the quasar luminosity 
function between z ~ 2 and z ~ 0.5 is not so strong in the 
model as it is in the dat a (here we have ma de the compari- 
son with the 2dF data bv lCroom et al.l2o'o3). This finding is 
not surprising because the figures in lKauffinaim fc Haehneld 



J2000h and lCattaneol i200ll) had already shown this limita- 
tion of the simplest scenario. The Durham group too has 
encountered the same problem when they have started try- 
ing and incorporating AGN into their semi-analytic model 
of galaxy formation (R. Malbon, private communication). 

4.2 The reference model 

The method to obtain a model that by construction repro- 
duces the strong evolution of the quasar population is to 
identify the parameters that are significantly different in 
high and low redshift AGN and to force a dependence of 
the black hole accretion rate on these parameters. We find 
that the main difference between z ~ 2 and z ~ 0.5 AGN 
is not in the gas fraction or in the potential well, but in the 
density of the gas in the central starburst that fuel the AGN. 
Therefore, we choose to explore a model in which the black 
hole accretion rate is 

M. oc p«^^^,M.burst. (8) 

Cattaneo, Haehnelt & Rees (1999) used a similar approach 
when they assumed that AM, oc (1 + 2)''AA/«burst with 
r; = 2 in order to obtain a reasonable agreement with the 
cosmic evolution of the comoving mass density of super- 
massive black holes. The difference is that here z does not 
appear explicitly. We find a reasonable agreement with the 
local black hole masses and with the luminosity function of 
quasars for a dependence with a power of = 0.5, a radia- 
tive efficiency of erad = 0.1 and model II for the fraction 
of obscured AGN (Figs. 5-6, model C, called the reference 
model hereafter). 

4.2.1 Black hole masses 

The blue (late-type galaxies) and red (early-type galaxies) 
point clouds in Fig. 6 show the results of the reference model 
at different redshifts while the squar es and the tria ngles with 
erro r bars are the mass es timates bv lMarconi fc H unt ( 20(^ 
and lHaring fc Rb3 (|20o3), and are the same in all four pan- 
els. In fact, they are the same as in Fig. 4. The black hole 
mass distribution at 2: = 2 is bimodal because two processes 
fuel AGN: mergers (predominant in elliptical galaxies) and 
bar instabilities (predominant in spiral and lenticular galax- 
ies). Low redshift starbursts are less dense and form objects 
with lower values of M, /Mbuigo. This mechanism generates 
scatter in the M, - Mbuigo relation. At low redshift the scat- 
ter becomes so large that it erases any trace of the original 
bimodality and creates a continuity between the products of 
bars and those of mergers. 

The Pburst factor not only introduces scatter, but also 
increases the slope of the Log(Af,) - Log(A/buigG) relation, 
since the densest starburst are also the most massive ones. 
Perhaps it is not coincidental the data contain a tilt in the 
same sense (the lines in Fig. 6 show the M, oc Mbuige relation 
found bvi Haring fc R ix 2004). 

Since iFerrarese fc Merrltl J2000l) and iGebhardt et alJ 

i2000l) . there has been considerable interest in the M. — a 
relation, where a is the velocity dispersion of the host bulge. 
This interest has been due to the discovery that the A/, — a 
relation contains less scatter than the M, — Lbui relation, 
where Lb ui is the optical luminos ity o f the host bulge. Nev - 
ertheless. iMarconi fc Hu'n3^l2003^ and lHaring fc Ri^ t004\ 
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F igure 5. The outpu ts of AGNICS (solid lines) for the six models listed in Table 1 compared with the 2dF quasar luminosity function 
of lCroom et alJ i2004) (points with error bars). The shaded areas show the error bars in the simulated luminosity functions. 
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Figure 6. The small points show the simulated black hole mass - bulge mass relation in a model in which M, = 
urst/kpc) "'^'^-^^burst (the reference model). Here Afcoldi ^burst and -A^*burst '^^^ th& gas mass, the Hernquist 
radius and the star formation rate of the starburst component. The blue points are for black holes and bulges of late-type galaxies. The 
red points are for black holes a nd bulges of early-ty pe galaxies. The open squares are the mass estimates bv lMarconi &: HuntI j20o3) 
and the filled triang les those bv lHaring fc Ri^3 i20n4h . The dashed-dotted lines corresponds to a 1.12 slope (the best fit to the data by 
iHaring fc Rb3l2004l) and the dashed lines show a scatter of ±0.3 dex around this relation. These lines are the same in all four panels. 
They refer to data at z = but they have also been shown in the diagrams at 2 > to ease the comparison. The results of the AGNICS 
model are shown for four different redshifts. The formal resolution of the simulation is 2 X IO^^A/q for the galaxy and 4 X IO'^A/q for 
the black hole. 



have shown that the M, — Mbui is as tight as the M. — u 
relation. Therefore, there is no reason to believe that the lat- 
ter is more fundamental than the former. However, we have 
chosen to show our simulated Mbuige — o" relation (Fig. 8). 
This appears to have a sh allower slope than that inf erred 
from the mass estimat es of iFerrarese fc Merritd i2000l) and 
iTremaine et alJ i2002h . but the difference is entirely at- 



tributable to the slope in the Mbuige — o" relation, which is 
shallower in GallCS than in the data (see the dashed line in 
Fig. 8). The result is not so bad if we consider that GallCS 
calculate radii, and therefore velocity dispersions, of bulges 
assuming that mergers conserve the mass and the total en- 
ergy, while neglecting the loss of mass and energy in tidal 
tails and ignoring any angular momentum consideration. 
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Figure 7. Cosmic evolution of the com oving number den sity of AGN at optical and X-ray wave l ength s. The points with error bars 
are the data. The optical data are from IWolf et all i2n0,'jl . The X-ray data are from lUeda et all i2003h . The lines are the results of 
our reference model. In the X-ray plot, we have two types of lines: the dashed lines corresponds to a bolometric correct ion of 32 as in 
thelElvis et al.l il994) median SED, the solid lines are computed with the luminosity-dependent bolometric correction in lMarconi et all 
J2004) . 



4-2.2 Evolution of the AGN population 

The fiducial model assumes that M, oc p{J,f^j,tMtburst- The 
scatter introduced by the Pbu^st factor is large enough that 
it allows to fit the luminosity function of quasars with trad = 
0.1 (model C of Fig. 5) without violating the constraints on 
black hole masses (Fig. 6). In Fig. 7, we take exactly the 
same model a nd compare its pr edictions with the data from 
COMBO-17 ("W olf et alj|2003ll . which probe fainter AGN 
and higher redshifts (left panel). In the right panel, we show 
the predictions of our model for the comoving number den- 
sity of AGN with 10*^ erg s"^ < L2-iokov < 3 x lO"*" erg s"^ 
and 3 x 10** erg s"^ < L2_iokcV < 10**ergs~\ The solid 
lines and the dashed lines are t he predictions d erived with 
the bolometric c orrections in iMarconi et alj ll2004ll and 
lElvis et al] (^^), res pectively. We com pare these predic- 
tions with the data bv lUeda et al.l (|20o3) (points with error 
bars). T he model that uses th e bolometric corrections pro- 
posed bv lMarconi et aO (|22q3) is in better agreement with 
the data. 

The fraction of obscured AGN in model C is computed 
with extinction law II, given by the dashed line in Fig. 3. 
We should clearly say that this law has been deliberately 
constructed to o btain the best possi ble agreement between 
model C and the'Croo m et alj (|2004|) data. A posteriori, the 
[Bargcr ct al. (2005) data have shown that this model is not 
inconsistent with our present knowledge of the type 2 frac- 
tion as a fu nction of l uminosity, since the d ata points from 
lUeda et al.l ^003) and lSzokolv et alj J2004) in Fig. 3 do not 
include the contribution of Compton-thick AGN to the type 
2 fraction (model D shows an alternative version of the ref- 
erence model, in which we only include the obscuration that 
is known to exist from AGN that are detected in X-rays, but 



do not have UV excess or broad line emission). One can be 
worried because the magnitude intervals on which the model 
and the data are compared at z ~ 2 and z ~ 1 only overlap 
at one point {Mb = —24.25). Reassuringly, at that point 
the agreement is good at both redshifts, but the comparison 
with the COMBO-17 data in Fig. 7 gives the best evidence 
that the success of our model is not a mere consequence of 
the change in the extinction law at the magnitude separat- 
ing high and low redshift data. Fig. 7 also shows that our 
model can reproduce a good agreement with the X-ray data. 
This agreement is independent of the obscuration model and 
supports the conclusion that the discrepancy between the 
blue luminosity function without obscuration and the opti- 
cal data is due to extinction. 

We introduced a steep increase in the type 2 fraction at 
low l uminosities to suppress the num ber of AGN with Mb > 
-24. lKauffmann fc Hachnclt ( 2000*) and lEnoki et al.l (l2003h 
dealt with the same problem by introducing a low circular 
velocity cut-off fEgs. I^Ht . Model F shows our results when 
we apply this solution in combination with the e xtinction 
law I derived from the data in iUoda ct al.^ i2003ll . The fit 
to the data is good at z ~ 2 and z ~ 0.5. The discrepancy 
at z ~ 1 occurs because all our AGN z ~ 1 are in large 
galaxy groups, where the low circular velocity cut-off has 
limited consequences, but this may be an effect of small 
number statistics due to cosmic variance together with the 
finite size of the computational box. The model with the low 
circular velocity cut-off gives a steeper slope in the Log(M,)- 
Log(Mbuige) relation than the model without this cut-off. 

Model E is an alternative version of model C in which 
we have dropped the Lboi < ^Edd condition. The excess of 
blue light from super-Eddington AGN is compensated by 
reducing the radiative efficiency from e^ad ~ 0.1 to erad = 
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Figure 8. Black hole and bulge mass plotted against the ve- 
locity dispersion of the host bulge at the half mass radius. The 
squ ares with error bars are the black hole mass estimates of 



iTremain ^^^ay j200j) while the d iamonds with error bars are 
those o nF^^ires^^ ^MerritjJjgQjjnh^I^ continuous solid line is 
a fit to the data"of l^^maineetaTl i |20o3) and has a slope of 4.02. 
The dashed line has a slope of 2.4. It is not a fit and is only shown 
to guide the eye. 



0.08. Model E fits the data even better that model C does, 
but we have chosen to present model C as our reference 
model because its assumptions are more standard. 



4-2.3 AGN activity and black hole mass 

Fig. 9 shows relation between the mass of the black hole and 
the luminosity of the AGN in model E (the reference model 
without the Lboi < ^Edd condition). This figure prompts 
two considerations. The first one is that, even when we re- 
move the Eddington limit, it is very difficult to find quasars 
brighter than 3— 4 Lsdd- The second one is the limited statis- 
tics deriving from the size of the computational box. At 
z = 2, there are 22 quasars with AIb < —24 and only 4 with 
Mb < -26. 

Fig. 10 presents the cosmic evolution of the mass func- 
tion of supermassive black holes and the fraction of black 
holes that are active at an Als < —22 level as a func- 
tion of black hole mass. These plots are identical for models 
C and E. The black hole mass function shows that black 
holes are still growing at 1 < 2; < 2, but the growth is 
small at 2: < 1. The evolution in the interval 1 < 2 < 2 is 
stronger at M. ~ 3 x 10* M© than it is at M. > W^Mq. 
The fraction of active black holes decreases by almost two 
orders of magnitude between z ~ 2 and 2 ~ 0. The most 
active black holes are objects of ~ 6 x 10* at 2 ~ 2 
and ~ 2 X 10* at 2 ~ 0.5. The two panels combined 
suggest a picture in which the most massive black holes 
form earlier while less massive black holes continue grow- 
ing to lower redshifts. This is consistent with the finding 
that the most massive galaxies form at higher redshifts 
('Monaco et al."2000'; "Granato e t 811120011: ICorbin fc Vaccal 
[2002 ; Cattaneo fc BernardL20 03i^ also known as down-sizing 
of galaxy formation or anti-hierarchical evolution of the 
baryons with respect to the dark matter, although Figs. 1 



and 10 show that hierarchical models of galaxy formation 
can in fact reproduce this behaviour. 



5 DISCUSSION AND CONCLUSION 

Observational and theoretical arguments (see the Introduc- 
tion) suggest that the same physical mechanisms are respon- 
sible for the growth of bulges and supermassive black holes. 
The main questions are what are these mechanisms and if 
they work one-way only (from the galaxy to the black hole) 
or both ways (through AGN feedback). 

Mergers and disc instabilities provide a path for trans- 
forming late-type galaxies into early-type galaxies while 
driving a sudden fuel supply into the galactic nucleus. The 
importance of this process is demonstrated by hydrody- 
namic simulations and observed in the real Universe. Semi- 
analytic models of galaxy formation have incorporated the 
merger model since the early 1990s and have achieved a sub- 
stantial degree of success (e.g. Fig. 1), but some discrepan- 
cies (overcooling at the centre of massive haloes, late forma- 
tion of massive galaxies, low SCUBA counts in relation to 
the blue light that comes out of local galaxies) remain. Crit- 
ics of semi-analytic models have argued that the great com- 
plexity of the models and the inherent large number of free 
parameters allow for the possibility that incorrect assump- 
tions may be reconciled with the data. After all, the Ptole- 
maic model allowed to compute the ephemerides with rea- 
sonable accuracy, but its basic assumptions were false. This 
is a real danger when models are developed to reproduce a 
small number of observational constraints. However, the in- 
creasing volume of astronomical data from sub-millimetre, 
infrared, optical and X-ray bands heavily outnumbers the 
free parameters in the hand of the simulators. The pres- 
ence of the same discrepancies in models developed indepen- 
dently and the difficulty that the modellers are encountering 
in solving these problems prove that the semi-analytic ap- 
proach is robust (although some relevant physics are still 
missing) . 

In this paper, we have presented AGNICS (Active 
Galactic Nuclei In Cosmological Simulations), a hybrid ap- 
proach that incorporates large cosmological N-body simu- 
lations, a semi-analytic model of galaxy formation and a 
scheme for the growth of supermassive black holes. We have 
used this model to investigate a scenario where mergers and 
disc instabilities drive massive gas inflows from galaxy discs 
into compact central starbursts, while a small fraction of 
this gas fuels the growth of a supermassive black hole. If 
Af, and Mtburst are simply related by a constant of propor- 
tionality, that leaves us with a tight Af, cx Af^^j^^, relation 
with (/ < 1. In a standard semi-analytic model, this simple 
scenario is not consistent with the strong evolution of the 
quasar population a.t z <. 2 observ ed in optical studie s of 
the quasar luminosity function (e.g. ICroom et aljl20o3 and 
references the rein). 

Recently, iBarger et al.l (|200^ have plotted the evolu- 
tion with redshift of the rest-frame 2-8 keV comoving en- 
ergy density production rate for AGN with L2_8kGv > 
10'*^ergs~^ and shown that this has a clear peak at 2 ~ 1. 
This result goes in the sam e sense of an earlier study by 
ICattaneo fc Bernard! feOOSTl . where we assumed that super- 
massive black holes form at the same time as the stellar 
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Figure 9. The black hole mass - blue luminosity relation in model E (the reference model without the Eddington limit) at four different 
redshifts. Each diagram is divided in six intervals of black hole mass. The row of numbers in the lower part of each diagram gives 
the number of black holes in each mass interval for the (150 Mpc)^ box with 256'' particles in the dark matter simulation. The formal 
resolution for the black hole mass is ~ 4 X V)' Mq. 



populations of their host galaxies and used this hypothesis 
to calculate the evolution with redshift of the total mass 
accreted by supermassive black holes per comoving volume 
and unit time. We found a peak at z ~ 1.6, but this can be 
an overestimate if some stars have formed before the black 
hole. 

The discrepancy between X-ray and optical data is 
attributed to the presence of an obscured AGN popula- 
tion (t^r£e_2_A^N2jwhM^^ at low luminosi- 
ties lBarger et all2005l: ISzokolv et alJ2004 lUeda et al.l2003l 
Low redshift AGN are less powerful and therefore more heav- 



ily obscured, with the consequence of amplifying the per- 
ceived evolution at optical wavelengt hs. 

Following th e met ho d of ISoltad (ES) and 
IChokshi fc Turneil lll992ll . lYu fc Tremaind 12002) esti- 



mated the mass per unit cosmic volume accreted by optical 
quasars over the life of the Universe, p.s, by using the 
equation 



(9) 



where (j}{LB ,t)dLB is the number density of quasars with 
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Figure 10. Left panel : cosmic evolution of the black hole mass function in the reference model. Right panel: the fraction of black holes 
with Mb < —22 as a function of black hole mass and redshift. In both panels, the solid, dotted, short-dashed and long-dashed lines 
correspond to z =2, 1, 0.5, 0, respectively. The formal resolution for the black hole mass is ~ 4 X 10^ Mq. 



blue luminosity between Lb and Lb + dLs at a time t af- 
ter the Big Bang and Cb = iboi/is (this equation uses 
Eg. El in the opposite d irection). For Srad = . 1, Cb = 11-8 
(as in lElvis et ai]ll994) and the iBovle et al.1 ll200(]ll model 
for the luminosity function of quasars, in which (j){LB,t) = 
4>*{l^ + l^), where I = LB/LB(t) a nd a and b are pa, - 
rameters determined from the data, IYu fc Tremaind (l2002l) 
found p,B — 2.1 X IO^MqMpc"^. They compared this value 
to their estimate of the local mass density of supermas- 
sive black holes, p. ~ (2.5 ± 0.4) x lO'^MoMpc"^ and con- 
cluded that the need for non-luminous accretion is lim- 
ited. There are two objections to this argument. The first 
one is that inserting <fi{LB,t) = 0*(/" -I- l^) into Eq. Q 
extrapolates the powe r law behaviour ou tside the mea- 
sured magnitude range. iBarger et al.l l^200!^^ argue that this 
leads to overestima ting p,b by a factor of ^ 1.8. Secondly, 
most other studies fSalucci et al.|[l999l : iMerritt fc Ferraresd 
120.01; Marconi ct al. 2004) favour values closer to p. ~ 4 x 
IO'^MqMpc"^. It seems plausible to conclude that 50 - 75% 
of bla ck hole accretion is optically obscured. ISazonov et alJ 
l)2004h analysed the combined infra-red and X-ray back- 
ground and reached a similar conclusion. 

With these considerations in mind, we ask ourselves if 
obscuration, which is more effective on the low redshift lumi- 
nosity function, can reconcile the simple mode M, oc M«burst 
with the data. The answer is that it cannot because model 
A in Fig. 5 does not fit the data even if it includes a large 
amount of obscuration (which may be an overestimate of the 
type 2 fraction at low luminosities). 

There are two possible solutions. The first one is 
that this problem is related to the other shortcomings of 
semi-analytic models of galaxy formation. Overcooling in 
massive haloes, late formation of elliptical galaxies, low 
number of high redshift sub-millimetre sources, delayed 



quasar epoch, all point to a scenario in which the most 
massive baryonic structures form too late. In recent times, 
there has been great interest on mechanical (jets, winds) 
and radiative (Compton heating) AGN feedback as a 
possi ble answer to this p h ysical p roblem (Tabor fc B innevI 
1991 iBinnevfcTaboi] fTooi iTucker fc DavidI |1997: 
Ciotti fc Ostriked Il997l I2001I: lOuilis et alJ 2001 : 
Revnolds et alJ 1200 ll 120021: iBasson fc Alexanded 12005 : 
'Om ma et alJ 120041: lOmma fc BinnevI 120041) . In this case, 
quasars and galaxy formation cease to exist as separate 
problems because we cannot solve one without solving the 
other. We plan to present and discuss such a model in a 
forthcoming publication of the AGNICS series. 

The second possibility, which we have studied in this 
paper, is to introduce a second parameter that breaks the 
strict proportionality between black hole accretion and star 
formation. This solution introduces a scatter and a tilt in 
the M, - Mbuige relation. 

A model in which the black hole accretion rate is en- 
hanced by a factor proportional to the square root of the gas 
density can reproduce the evolution of the quas ar population 
over the entire redshi ft range covered b y 2dF iCroom et alJ 
12004") and Combo- 17 ("Wolf et al.''200^ data both at optical 
and X-ray l)Ueda et al...2003,) wavelengths, without need to 
deviate from the standard radiative efficiency of erad = 0.1 
and in a manner that is consistent with the observational 
constraints on the normalisation, the slope and the scatter 
of the M. - Afbuigc relation. This model produces a final 
black hole mass density of p, ~ 4.4 x IO^MqMpc"^, very 
close to the most rec ent estimate of 4.5 x IO^MqMpc"^ by 
iMarconi et al.l i2004h (see the discussion above). Black holes 
with M. ^ 1 — 2 X 10* M© are in bulge dominated galaxies 
and have mostly grown through mergers. At M, < 10* 
disc instabilities are an important growth path. The good 
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Figure 11. The point cloud shows the relation between the bolo- 
metric luminosity of the AGN and the bolometric luminosity of 
the starburst component due to star formation only. This rela- 
tion has been plotted using all the AGN present in the refer- 
ence simulation at all timesteps. The solid line corresponds to 

-^'AGN = -^'starburst- 



and sub-millimetre properties of AGN is left to a future pub- 
lication of the AGNICS series. We can anticipate that we 
have checked the most fundamental constraints. E.g. even if 
we make the extreme assumption that all the AGN power 
in the Universe is absorbed and reemitted in the far in- 
frared by dust in the host galaxies, that could account for 
< 1/3 of the number of SCUBA sources counted at a flux 
of Sssoum > 2mJy. The second prediction of our model is 
that, for a given bulge mass, the most massive black holes 
are found in the bulges with the old est stellar populations. 
Observations iMorrifield et al.ll200Cl) seem to confirm this 
second prediction, too, but it is important to have more data 
to make this conclusion more secure. Moreover, low redshift 
mergers have much lower gas fraction than high redshift 
mergers. That can be another explanation why bulges with 
younger stellar populations have a lower black hole mass 
fraction. The first prediction is thus a cleaner test of our 
assumption. 
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agreement with the X-ray data support the conclusion that 
the high type 2 fraction that we need to invoke to fit the 
optical data at low luminosities has a physical basis. Allow- 
ing or preventing radiation at Lboi > ^Edd does not any 
significant difference, which cannot be compensated by a 
small change of the value of the radiative efficiency trad. It 
is difficult to say to which extent this model is physical and 
to which extent the factor oc pburst simply compensates the 
shortcomings mentioned above. Up to a point, it certainly 
does, but it would be a coincidence if the dependence that 
compensates a shortcoming of the galaxy formation model 
and gives the appropriate evolution of the quasar luminosity 
function was also the dependence that introduces the right 
scatter in the black hole mass distribution (and tilts the dis- 
tribution of an amount that does not exceed the constraints 
on the power of the M, - Mbuigc relation). 

The solution proposed in this paper, which breaks the 
exact proportionality between black accretion and star for- 
mation in favour a higher accretion rate in dense high red- 
shift starbursts, makes two observable predictions. The first 
one is that the quasar contribution to the total bolometric 
luminosity is stronger in the most powerful, denser star- 
bursts (Fig. 11). This is a prediction that can be tested 
by inspecting the SEDs of ultra-luminous infra-red galax- 
ies (ULIRGs) for AGN features and by studying how the 
presence and the strength of these features in crease with 
the to tal bolometric luminosity of these objects. iTran et alJ 
i200ll) studied the mid-infrared spectra of 16 ULIRGs and 
found a transition from mostly starburst-powered to mostly 
AGN-powered objects at a luminosity of ~ 2 — 3 x 10^^ Lq 
in reasonable agreement with the prediction of the reference 
model (Fig. 11). The systematic discussion of the infra-red 
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